*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* Program for creating data and specifying model estimation		~~~~
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

capture program drop sim
program define sim, rclass
	version 14.1
	syntax, obs(numlist) waves(numlist)
	clear

*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* Creating the data												~~~~
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

* First create data structure
mat C =(1, $corrXY,$corrXZ \ $corrXY,1,$corrYZ \ $corrXZ,$corrYZ,1)
corr2data X Y Z, n($obs) corr(C)

* Then expand data to number of waves as specified in master file	
gen id = _n
expand = $waves
bysort id: generate wave = _n
replace wave = wave-1
foreach var of varlist X Y Z {		// missing values for t>1
	replace `var' = . if wave >0	
}
xtset id wave
reshape wide X Y Z , i(id) j(wave)

* Then generate data for t1		
gen x0 = X0 + rnormal()				// xt1 + error term
gen y0 = Y0 + rnormal()				// yt1 + error term
gen z0 = Z0 + rnormal()

* Then run flexible data-generating process for further waves
foreach num of numlist 1/$waves{			
gen x`num' = $autox*x0 + rnormal() + $endo* y`=`num'-1' + ($unob*z0) 
gen y`num' = $autoy*y0 + rnormal() + 			///
((1-$lambdax)*$betax*(x`num')) + (				///
$lambdax*$betaxlag*(x`=`num'-1')) + ($unob*z0)
}	
			
reshape long x y ylone ylhalf ylzero,i(id) j(wave)
keep if wave >0		

*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
* Estimating models												~~~~
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

* Pooled OLS (POLS)
local model = "POLS"					
reg y x, vce(robust)					// POLS: contemp. effect	
	return scalar `model'con			= _b[x]
	return scalar `model'conse 			= _se[x]
	return scalar `model'conlb95 		= _b[x] + 1.96 * _se[x]
	return scalar `model'conub95 		= _b[x] - 1.96 * _se[x]
reg y L.x, vce(robust)					// POLS: lagged effect
	return scalar `model'lag 			= _b[L.x]
	return scalar `model'lagse 			= _se[L.x]
	return scalar `model'laglb95 		= _b[L.x] + 1.96 * _se[L.x]
	return scalar `model'lagub95 		= _b[L.x] - 1.96 * _se[L.x]
reg y x L.x, vce(robust)				// POLS: both effects
	return scalar `model'both_con 		= _b[x]
	return scalar `model'both_conse 	= _se[x]	
	return scalar `model'both_conlb95 	= _b[x] + 1.96 * _se[x]
	return scalar `model'both_conub95 	= _b[x] - 1.96 * _se[x]		
	return scalar `model'both_lag 		= _b[L.x]
	return scalar `model'both_lagse 	= _se[L.x]			
	return scalar `model'both_laglb95 	= _b[L.x] + 1.96 * _se[L.x]
	return scalar `model'both_lagub95 	= _b[L.x] - 1.96 * _se[L.x]
	
* Random Effects (RE)
local model = "RE" 						
xtreg y x, re vce(robust)				// RE: contemp. effect
	return scalar `model'con			= _b[x]
	return scalar `model'conse 			= _se[x]
	return scalar `model'conlb95 		= _b[x] + 1.96 * _se[x]
	return scalar `model'conub95 		= _b[x] - 1.96 * _se[x]
xtreg y L.x, re vce(robust)				// RE: lagged effect
	return scalar `model'lag			= _b[L.x]
	return scalar `model'lagse 			= _se[L.x]
	return scalar `model'laglb95 		= _b[L.x] + 1.96 * _se[L.x]
	return scalar `model'lagub95 		= _b[L.x] - 1.96 * _se[L.x]		
xtreg y x L.x, re vce(robust)			// RE: both effects
	return scalar `model'both_con 		= _b[x]
	return scalar `model'both_conse 	= _se[x]	
	return scalar `model'both_conlb95 	= _b[x] + 1.96 * _se[x]
	return scalar `model'both_conub95 	= _b[x] - 1.96 * _se[x]		
	return scalar `model'both_lag 		= _b[L.x]
	return scalar `model'both_lagse 	= _se[L.x]			
	return scalar `model'both_laglb95 	= _b[L.x] + 1.96 * _se[L.x]
	return scalar `model'both_lagub95 	= _b[L.x] - 1.96 * _se[L.x]	
	
* Fixed Effects (FE)	
local model = "FE" 						
xtreg y x, fe vce(robust)				// FE: contemp. effect
	return scalar `model'con			= _b[x]
	return scalar `model'conse 			= _se[x]
	return scalar `model'conlb95 		= _b[x] + 1.96 * _se[x]
	return scalar `model'conub95 		= _b[x] - 1.96 * _se[x]
xtreg y L.x, fe vce(robust)				// FE: lagged effect
	return scalar `model'lag			= _b[L.x]
	return scalar `model'lagse 			= _se[L.x]
	return scalar `model'laglb95 		= _b[L.x] + 1.96 * _se[L.x]
	return scalar `model'lagub95 		= _b[L.x] - 1.96 * _se[L.x]		
xtreg y x L.x, fe vce(robust)			// FE: both effects
	return scalar `model'both_con 		= _b[x]
	return scalar `model'both_conse 	= _se[x]	
	return scalar `model'both_conlb95 	= _b[x] + 1.96 * _se[x]
	return scalar `model'both_conub95 	= _b[x] - 1.96 * _se[x]		
	return scalar `model'both_lag 		= _b[L.x]
	return scalar `model'both_lagse 	= _se[L.x]			
	return scalar `model'both_laglb95 	= _b[L.x] + 1.96 * _se[L.x]
	return scalar `model'both_lagub95 	= _b[L.x] - 1.96 * _se[L.x]	
	
* First-Differences (FD)	
local model = "FD" 						
reg D.y D.x, vce(robust)				// FD: contemp. effect
	return scalar `model'con			= _b[D.x]
	return scalar `model'conse 			= _se[D.x]
	return scalar `model'conlb95 		= _b[D.x] + 1.96 * _se[D.x]
	return scalar `model'conub95 		= _b[D.x] - 1.96 * _se[D.x]
reg D.y L.D.x, vce(robust)				// FD: lagged effect
	return scalar `model'lag			= _b[L.D.x]
	return scalar `model'lagse 			= _se[L.D.x]
	return scalar `model'laglb95 		= _b[L.D.x] + 1.96 * _se[L.D.x]
	return scalar `model'lagub95 		= _b[L.D.x] - 1.96 * _se[L.D.x]		
reg D.y D.x L.D.x, vce(robust)			// FD: both effects
	return scalar `model'both_con 		= _b[D.x]
	return scalar `model'both_conse 	= _se[D.x]	
	return scalar `model'both_conlb95 	= _b[D.x] + 1.96 * _se[D.x]
	return scalar `model'both_conub95 	= _b[D.x] - 1.96 * _se[D.x]		
	return scalar `model'both_lag 		= _b[L.D.x]
	return scalar `model'both_lagse 	= _se[L.D.x]			
	return scalar `model'both_laglb95 	= _b[L.D.x] + 1.96 * _se[L.D.x]
	return scalar `model'both_lagub95 	= _b[L.D.x] - 1.96 * _se[L.D.x]				

* Arellano-Bond (AB)	
local model = "AB" 						
xtabond2 y L.y  x, gmm(L.y x) noleveleq robust  		// Arellano-Bond: contemp. effect
	return scalar `model'con 			= _b[x]
	return scalar `model'conse 			= _se[x]
	return scalar `model'conlb95 		= _b[x] + 1.96 * _se[x]
	return scalar `model'conub95 		= _b[x] - 1.96 * _se[x]			
xtabond2 y L.y  L(1/1).x, gmm(L.y L.x) noleveleq robust // Arellano-Bond: lagged effect
	return scalar `model'lag 			= _b[L.x]
	return scalar `model'lagse 			= _se[L.x]	
	return scalar `model'laglb95 		= _b[L.x] + 1.96 * _se[L.x]
	return scalar `model'lagub95 		= _b[L.x] - 1.96 * _se[L.x]
xtabond2 y L.y  L(0/1).x, gmm(L.y L.x) noleveleq robust 
	return scalar `model'both_con 		= _b[x]
	return scalar `model'both_conse 	= _se[x]	
	return scalar `model'both_conlb95 	= _b[x] + 1.96 * _se[x]
	return scalar `model'both_conub95 	= _b[x] - 1.96 * _se[x]		
	return scalar `model'both_lag 		= _b[L.x]
	return scalar `model'both_lagse 	= _se[L.x]			
	return scalar `model'both_laglb95 	= _b[L.x] + 1.96 * _se[L.x]
	return scalar `model'both_lagub95 	= _b[L.x] - 1.96 * _se[L.x]		

* Cross-Lagged SEM with Fixed Effects (ML-SEM)	
local model = "MLSEM" 			
capture xtdpdml y, pre(x) vce(robust)		// ML-SEM: contemp. effect
	if e(converged)== 1 {					// Only converged models
		return scalar `model'con 			= _b[y3:x3]
		return scalar `model'conse 			= _se[y3:x3]
		return scalar `model'conlb95 		= _b[y3:x3] + 1.96 * _se[y3:x3]
		return scalar `model'conub95 		= _b[y3:x3] - 1.96 * _se[y3:x3]		
	}				
capture xtdpdml y, pre(L.x) vce(robust)		// ML-SEM: lagged effect
	if e(converged)== 1 {					// Only converged models
		return scalar `model'lag 			= _b[y3:x2]
		return scalar `model'lagse 			= _se[y3:x2]
		return scalar `model'laglb95 		= _b[y3:x2] + 1.96 * _se[y3:x2]
		return scalar `model'lagub95 		= _b[y3:x2] - 1.96 * _se[y3:x2]	
	}
capture xtdpdml y, pre(x L.x) vce(robust)	// ML-SEM: both effects
	if e(converged)== 1 {					// Only converged models
		return scalar `model'both_con 		= _b[y3:x3]
		return scalar `model'both_conse		= _se[y3:x3]
		return scalar `model'both_conlb95 	= _b[y3:x3] + 1.96 * _se[y3:x3]
		return scalar `model'both_conub95 	= _b[y3:x3] - 1.96 * _se[y3:x3]					
		return scalar `model'both_lag		= _b[y3:x2]
		return scalar `model'both_lagse		= _se[y3:x2]
		return scalar `model'both_laglb95 	= _b[y3:x2] + 1.96 * _se[y3:x2]
		return scalar `model'both_lagub95 	= _b[y3:x2] - 1.96 * _se[y3:x2]	
	}	
end

